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We consider a system of two discrete nonlinear Sciirodinger equations, coupled by nonlinear and 
linear terms. For various physically relevant cases, we derive a modulational instability criterion for 
plane-wave solutions. We also find and examine domain-wall solutions in the model with the linear 
coupling. 

I. INTRODUCTION 

Modulational instabilities (Mis) have a time-honored history in nonlinear wave equations. Their occurrences span 
areas ranging from fluid dynamics l] (where they are usually referred to as the Benjamin- Feir instability) and nonlinear 
optics (2, i3] to plasma physics ^ . 

While earlier manifestations of such instabilities were studied in continuum systems 0, , in the last decade the 
role of the MI in the dynamics of discrete systems has emerged. In particular, the MI was analyzed in the context of 
the discrete nonlinear Schrodinger equation , a ubiquitous nonlinear-lattice dynamical model 0, ll| • More recently, 
it was studied in the context of weakly interacting trapped Bose-Einstein condensates (BECs), where the analytical 
predictions based on the discrete model M. were found to be in agreement with the experiment [T^ (see also the recent 
works 0,0,0], and recent reviews in |l4lT5l|'). Additionally, the development of "discrete nonlinear optics" (based 
on nonlinear waveguide arrays) has recently provided the first experimental observation of the MI in the latter class 
of systems [l^ . 

An extension of these recent works, which is relevant to both BECs in the presence of an optical-lattice potential 
[itL IT^ |2iI| and nonlinear optics in photorefractive crystals l2^, is the case of multi-component discrete fields. 
These can correspond to a mixture of two different atomic species, or different spin states of the same atom, in BECs 
|23|| , or to light waves carried by different polarizations or different wavelengths in optical systems 24] . 

The corresponding two-component model is based on two coupled discrete nonlinear Schrodinger (DNLS) equations, 

= -rfl(wi,n+l + - 2ui,„) + (Sll]ui,„]^ + Sl2 ]u2,n ] ^)ui,n + CU2,n, 



dt 

.du2,n 



-d2(u2,n+l + W2,n-1 " 2^2, „) + (si2]ui,n]^ + S22 ]u2,n ] ^)"2,n + CUl,n- (1) 



This is a discrete analog of the well-known model describing nonlinear interactions of the above-mentioned light 
waves through self-phase-modulation and cross-phase-modulation (XPM) In optics, Eqs. (Q) describe an array 

of optical waveguides, the evolution variable being the propagation distance z (rather than time t). The choice of 
the nonlinear coefficients in the optical models is limited to the combinations sn = S22 = 3si2/2 for orthogonal 
linear polarizations, and sn = S22 = S12/2 for circular polarizations or different carrier wavelengths. In BECs, 
the coefficients Sjk in Eqs. are related to the three scattering lengths Ujk which account for collisions between 
atoms belonging to the same [ajj) or different (a^fe, j ^ k) species; in that case, > (a^- < 0) corresponds 
to the repulsive (attractive) interaction between the atoms. The linear coupling between the components, which is 
accounted for by the coefficient c in Eqs. is relevant in optics for a case of circular polarizations in an array of 
optical fibers with deformed (non-circular) cores, or for linear polarizations in an array of twisted fibers j24j . On the 
other hand, in the BECs context, c represents the Rabi frequency of transitions between two different spin states in 
a resonant microwave field |2(l. Wv{ . 

Stimulated by the experimental relevance of the MI in discrete coupled systems, the aim of the present work is 
to develop a systematic study of the instability in the two-component dynamical lattices. We give an analytical 
derivation of the MI criteria for the case of both the nonlinear and linear coupling between the components. As a 
result of the analysis, we also find a novel domain-wall (DW) stationary pattern in the case of the linear coupling. 
The presentation is structured as follows: In the following sections we derive plane- wave solutions and analyze their 
stability, corroborating it with a numerical analysis of the stability intervals. We do this for the model with the linear 
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coupling in section II and for the one with the purely nonlinear coupling in section III. In section IV, we examine DW 
states in the linearly coupled lattices. Finally, the results and findings are summarized in section V. 



II. THE MODEL WITH THE LINEAR COUPLING 



We look for plane-wave solutions in the form 

u jn ^ Aj ejip[i{qjn - ujjt)] , j = l,2. (2) 
The linear coupling imposes the restrictions qi = q2 = q and cji = lj2 = t^- Inserting Eq. |(5Jl into Eqs. |^ yields 

ujAi = -2di{cofiq-l)Ai + {snAl + si2Al)Ai+cA2, 

UJA2 = -2d2{cOSq-l)A2 + {si2Al+S22Al)A2+cAi. (3) 

In the particular, but physically relevant, symmetric case, with di = d2 and sn = S22, it follows from here that the 
amplitudes Ai^2 obey an equation [(sn — si2)Aiyl2 — c] {Af — A2) = 0, hence either of the following two relations 
must then be satisfied: 

Al = ±A2; (4) 



A1A2 = . (5) 

Sn - S12 



In the case of Eq. Q), a nonzero solution has Ai = zLy/{2di{cosq — 1) + ^ c) / (sn + S22). It exists with 
Sn + S12 > 0, provided 2(ii(cos(7 — 1) + w =F c > 0, and with sn + S12 < 0, if 2di(cosg — 1) + uj ^ c < 0. When 
Sll = —512, one obtains solutions of the form {Ai, A2) = {±A, A) with arbitrary A and — 2(ii(cos q — 1) — uj ^ c. 

On the other hand, from Eq. 0, one finds that 

2 uj + 2c?i(cosg — 1) 
Al = 




^Aicosq^V A^ (6) 

Sn / (sii - si2)^ 

under the restriction that this expression must be positive. When c 7^ 0, solutions of this type exist as long as 
[u> + 2di(cos(7 — 1)] Sn > (as the term under the square root is smaller in magnitude than the one outside) and the 
argument of the square root in © is non-negative. The first condition implies lu > — 2cii(cosg — 1), for sn > 0, and 
uj < —2di{cosq — 1), for sn < 0. In the case S12 — 0, Eq. (jSJ takes the simpler form 



A? = (2sii) ^ uj + 2di(cos(j - 1) ± + 2di(cos(? - 1))2 - 4c^ . (7) 

For both S12 = and S12 = 2sii, it is necessary to impose the condition \uj + 2cii(cosg — 1)1 > 2c for the solutions to 
be real. We note that these relations are similar to those derived in Ref. where linearly and nonlinearly coupled 
systems of continuum NLS equation were considered. 

It is also worth noting that the existence of two distinct uniform states with a fixed product from Eq. |(SJ) suggests a 
possibility of a domain- wall (DW) solution in the model with the linear coupling. DW solutions in nonlinearly coupled 
discrete nonlinear Schrodinger equations were examined in Ref. "25] (following an analogy with the continuum ones 
of Ref. H^). However, the present case is different, as both uniform states have non-vanishing amplitudes in both 
components, and, as seen from Eq. (jSJ, the DWs may exist only if c 0. This possibility is examined in more detail 
in section IV. 

To examine the stability of the plane waves, we substitute 

Ujn{x, t) = [Aj + Bjn{x, t)] exp[i(qn - uot)] (8) 

into Eqs. to obtain a system of two coupled linearized equations for the perturbations Bj{x,t). Furthermore, 
assuming a general solution of the above-mentioned system of the form 

Bjn = aj cos{Qn — fit) + i(3j sin(Qn — fit), (9) 
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where Q and fl are the wavenumber and frequency of perturbation, we arrive at a set of four homogeneous equations 
for ai, (3i, a2 and ^2- The latter have a nontrivial solution if Q and satisfy the dispersion relation 



(fl — 2di sin Q sinqf' — ( 2di cos q{cos Q — 1) + c-^ ) [ 2di cosq(cos Q — 1) + c—^ — 2siiA^ 



A, 



{VL — 2d2 sin Q sin g)"^ — ( 2^2 cosg(cos(3 — 1) + j \^d2 cos (/(cos Q — 1) + — 2S22A. 



— 2c(2si2 A2 + c) (il — 2c?i sin Q sin q) {fl — 2d2 sin Q sin q) - 
X ^2^2 cos 9(cos Q-l) + - 2s22Al \ - {2S12A1A2 + cf 



Ao 



2di cos g(cos Q — 1) + c— 2siiA'l 

Ai 



2di cosg(cos(5 — 1) + j ('2d2 cos(7(cos(5 — 1) + ] — 



= 0. (10) 



(11) 



Note that in the absence of coupling, i.e., c — S12 — 0, we obtain a known relation [g 

(il — 2dj sin Q sin qj)"^ — 2dj cos qj (cos Q — 1) (2c?j cos qj (cos Q — 1) — 2sjjj4^) , 

which gives the MI condition if the right-hand side becomes negative, i.e., 

dj cosgj(2dj cosq^ sin^ (Q/2) + SjjA^j) < 0. (12) 

We will now consider the particular case di ~ d2 = d, the general case being technically tractable but too involved. 
Then, the dispersion relation (|10|) takes the from 

{n - 2d sing sin q)4 - {Ki + K2 + K3){n - 2d sin g sin g)^ + K1K2 - = 0, (13) 



where 



2dcosg(cosg — 1) + c— ^ ) ( 2dcosg(cosg — 1) + c— ^ — 2siiAl 
^1/ V Ai 



K2 = ^2d cos g(cos g - 1) + ^2d cos (?(cosg - 1) + - 2522^2 

Ks. = 2c(2si2^iA2 +c). 



2dcosq(cosg — 1) 



c^-2s,,Al 



-{2si2AiA2 + cf 



2d cos (/(cos g — 1) 



2d cos q (cos g 
A2 



- 2522^2 

A2 



'A, 



2d cos (/(cos g — 1) 



'A2 



It immediately follows from Eq. 113|) that, to avoid the MI, both solutions for (fi — 2d sing sin (/)^ 2 should be positive. 
Taking into account the binomial nature of the equation, it is concluded that the spatially homogeneous solution is 
unstable if either the sum S = Ki + K2 + K-^ or the product 11 = K1K2 — K4 of the solutions is negative: 



K1 + K2 + K3 < 0; 
K1K2 -Ki < 0. 



(14) 
(15) 



To proceed further, we may fix the value of the perturbation wavenumber, to investigate in what parameter 
region it would give rise to the MI. We illustrate this approach in Figs. ^ and 12 in which we fix g = tt and 
sii = S22 = ^1 = j42 = di d2 = 1 and vary c and S12 (the coefficients of the linear and XPM coupling), to examine 
their effect on the stability interval. From Eq. H12|l we see that for these values of the parameters the modulational 
unstable region is tt/2 < q < 2tt/3 = 2.0945. It can be inferred from the figures that c may widen the MI interval by 
decreasing its lower edge. On the other hand, S12 has a more complex effect: while making the instability interval 
larger by increasing its upper edge (until it reaches tt) , it may also open MI bands within the initially modulationally 
stable region (see, e.g., the lower panel of Fig. O. 

For reasons of completeness, we also illustrate nonlinear development of the MI in the coupled nonlinear lattices in 
some typical examples. In particular, the role of the MI in generating large-amplitude excitations in the presence of 
the linear coupling only {si2 = 0) is illustrated in the left panel of Fig. |31 for c = 0.5. In the right panel of the figure, 
the MI in the presence of both linear and nonlinear coupling (c = 0.5, S12 = 2/3) is shown. It is readily observed that 
the nonlinear coupling enhances the instability growth rate and, hence, the MI sets in earlier. 




FIG. 1: The figure shows, for S12 — 0, the cases of c = (top panel; unstable for n/2 < q < 2.0945), c = 0.25 (middle panel; 
unstable for 1.4455 < q < 2.0945) and c = 0.5 (bottom panel; unstable for 1.318 < q < 2.0945). The solid line shows the sum 
E and the dashed line the product IT of the solutions of Eq. 11311 . The instability takes place in intervals of the wavenumber q 
of the unperturbed plane-wave solution where either E or 11 (or both) are negative. 
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FIG. 2: Same as Fig. 0but fixing c — 0.25 and varying S12. The figure shows the cases of S12 — 2/3 (top panel; unstable 
for 1.4455 < q < 2.556), S12 = 1 (middle panel; unstable for n/2 < q < n), and S12 — 2 (bottom panel; unstable for 

0. 8955 < g < 1.4455 and 7r/2 < g < tt). 

III. THE CASE OF THE PURELY NONLINEAR COUPLING 

A physically relevant case that we also wish to consider here is the one with only the nonlinear coupling present, 

1. e., c = 0. The corresponding dispersion relations read 

LUj = -2dj{coaqj - 1) + sjiAf + Sj2^i j 1, 2. (16) 
To study the stability of the plane waves in this case, we use Eq. JHJ as before, and obtain the following dispersion 
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FIG. 3: Left panels: Spatio-temporal (n, t) contour plots of the squared absolute value of one of the two components of the 
solution (top) and temporal evolution of the maximum of the solution's amplitude (bottom), c — 0.5, S12 — and q = n/2. 
Right panels; Same as the left ones but for c — 0.5, S12 — 2/3 and q — -k 12. 



relation for the perturbation wavenumber and frequency 

[(il — 2d\ sinQsingi)'^ — 2c?i cos(7i(cos(3 — 1) {2d\ cos(7i(cos(3 — 1) — 2siiylJ)] 
X [(il — 2^2 sin Q sin 52)^ — 2c?2 cos (72 (cos Q — 1) (2^2 cos 52 (cos Q — 1) — 2522^3)] 
- 4(2si2AiA2)^difi2CosgiCOsg2(cos(5- 1)^ = 0. (17) 

When S12 = 0, the known result 1)12(1 for the one-component case is easily retrieved. If we let di = d2 = d and 
qi — q2 = q, then Eq. ()17|l is simplified as follows: 

{n - 2d sin Q sing)* - 2/^5 (f7 - 2d sin Q sin g)^ + Kq = 0, (18) 

where 

— 2dcosq{cosQ — 1) (2dcosq{cosQ — 1) — {siiAl + S22A2)) , 
Ke = (2dcosg(cos(3 - 1))^ ((2dcos g(cos Q - 1))^ - 2{siiAl + S22Al)2dcosq{cosQ - 1) + 4:AlAl(siiS22 - s^)) 

(19) 

We will follow the lines of the analysis outlined above for the case of the linear coupling. Examining the product 
and the sum of the roots of Eq. (|18|) for (51 — 2(isin(3sing)^, which now read S = and 11 = Kq, we arrive at the 
following MI conditions: 

K5 < 0, or Ke < 0. 

The latter can be rewritten, respectively, as: 



-(siiAf + 522^2) < 2dcos(q)sin" < 0; (20) 

< -4dcos((j)sin2 ( ^ ) < (21) 



with K± = siiAf + S22A2 ± y/ (sii^i + 522^2)^ ^ 4Af ^2(511322 — 512)- It is readily observed, in this case as well, 
that the coupling between the two components tends to expand the band of the MI wavenumbers with respect to the 
single-component case. The effect of the variation of S12 in this case, for fixed su — S22 = Ai = A2 = di = d2 = 1 



6 



























4 


















- 


1—1 


2 


- 
















, - - 



















-^^^ 


























-2 








- 


- -r ' 












1.2 


1.4 


1.6 1. 


8 2 


2.2 


2.4 


2.6 


: \ 


3 




5 








































G 







\ 






































-5 






















1.2 


1.4 


1.6 1 


8 2 


2.2 


2.4 


2.6 


2.8 


3 









C 






















-10 




















w 


-20 
























1.2 


1.4 


1.6 1. 


8 2 


2.2 


2.4 


2.6 


2.8 


3 



q 



q 



FIG. 4: Same as Fig. |5]but setting c = and varying sn. The figure shows the cases of si2 = 2/3 (top panel; unstable for 
7r/2 < q < 2.5550), Si2 = 1 (middle panel; unstable for 7r/2 < q < n), and si2 = 2 (bottom panel; unstable for 1.0472 < g < vr). 



and Q = TT, is shown in Fig. 01 Notice also that Eqs. (|20|l and H21|l suggest multi-component generalizations of the 
MI criteria. 

When sii — S22 = s, Ai = A2 = A and s, S12 > 0, the MI conditions given by Eqs. ()20f) - ()21|l can be written in a 
compact form, 

- sA"^ < 2dcos{q) sin^ (Q/2) < 0; (22) 
(5-512)^2 < -2dcos((7)sin2(Q/2) < (s + si2)yl^ (23) 



From Eq. (|^ . we obtain that (for Q 



TT / sA 



12 



2 <q<n- arccos ( — ) , (24) 



while similarly from Eq. H23|) . it follows that 



i' {s-si2)A^ \ f {s + si2)A^ 
TT — arccos — 1 < q < n — arccos — I . (25) 



2d " \ 2d 



Examining the latter relations in more detail, we conclude that, for s > S12, the region of unstable wavenumbers is 

tt/2 < q < TT — arccos j while, for s < S12, Eq. (|25|l contains the interval of unstable q. This is illustrated 

for the case with A = ii = s = lasa function of S12 in Fig. O 

A typical example of the simulated development of the MI in the case of purely nonlinear coupling is shown in Fig. 
Elfor S12 = 2/3 and q = tt/2. 



IV. DOMAIN WALLS IN THE SYSTEM WITH THE LINEAR COUPLING 

The expressions (O for the amplitudes of the plane-wave solutions suggest a novel possibility: Focusing more 
specifically on the so-called anti-continuum limit of d = (we use di — d2 = d) and following the path of '25*1, we 
can construct domain-wall (DW) solutions that connect the homogeneous state given by Eq. Q with its "conjugate" 
state of A2 — c/{Ai(sii — S12)). Such a solution can then be continued for finite coupling, to examine its spatial 
profile and dynamical stability. 

An example of such a solution, for the case where the nonlinear coupling is absent, is given in Fig. It is observed 
that the solution is stable for d < 0.033, and it becomes unstable due to a cascade of oscillatory instabilities (through 
the corresponding eigenvalue quartets) for larger values of d. 
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FIG. 5: The threshold wavenumbers q for the modulational instabihty vs. si2, as found from Eqs. H24|l and H25^ (for the case 
of A = d = s = 1). 




FIG. 6: Same as Fig. |21 but for si2 — 2/3 and c — (purely nonlinear coupling); q = 7r/2. 



We have also examined the evolution of the DWs when they are unstable. A typical result is displayed in Fig. El It 
is seen from the bottom panel, which shows the time evolution of the solution's maximum amplitude, that growth of 
the oscillatory instability eventually destroys the configuration, through "lattice turbulence" . This apparently chaotic 
evolution can be attributed to the mixing of a large number of unstable eigenmodes. The dynamics remain extremely 
complex despite the eventual saturation of the instability. 

The increase of the nonlinear coupling constant Si2 reduces the stability window of these solutions. In particular, 
the case of si2 = 0.3 is shown in Fig. IHlin which, the stability window has shrunk to d < 0.012. 

One can also consider a modified DW, where an extra site between the two domains has equal or opposite amplitudes 
of the two fields. We have checked that such solutions are always unstable (due to the presence of real eigenvalue 
pairs), for all values of d. Still more unstable (with a larger number of unstable eigenvalues) are more sophisticated 
DW patterns, with additional intermediate sites inserted between the two domains. 
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FIG. 7: The left panel (top subplot) shows the evolution of the norm of the DW solution Uj„ = cxp{—iiot)vjn, as a function of 
the inter-sito coupling constant d using continuation from the anti- continuum limit, d — Q. The bottom subplot shows the most 
unstable eigenvalue of small perturbations around the solution as a function of d. The DW becomes unstable for d > 0.033. 
The profiles of the solution (top subplots) and the respective linear-stability eigenvalues (bottom subplots) are shown in the 
right panel for d = 0.01 (stable; left subplots) and d = 0.04 (unstable; right subplots). The parameters are si2 = 0, c = 0.25, 
s = u> = 1. 




FIG. 8: The top panel shows a grayscale image of the spatio-temporal evolution of the intensity of the first field, |win(i)l^- The 
middle panel shows the same for the second field. The bottom panel displays the time evolution of the amplitude (the spatial 
maximum) of the first field. 



V. CONCLUSIONS 



In this work, we have examined the extension of the modulational instability (MI) concept to the case of multiple- 
component discrete fields. We have shown, in a systematic way, how to formulate the linear MI equations and how 
to extract the MI criteria. We have also followed the dynamical evolution of the instability by means of direct 
simulations, and have identified the effects of the linear and nonlinear couplings on the range of modulationally 
unstable wavenumbers. In particular, we have demonstrated that the joint action of the two couplings may give rise 
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FIG. 9: Same as Fig. but for S12 = 0.3. The instability sets in for d > 0.012. 



to noteworthy features, such as opening of new MI bands on the wavenumber scale. 

Additionally, the identification of a pair of conjugate uniform solutions in the two-component model has prompted 
us to examine domain-wall (DW) solutions between such states. We were able to demonstrate that the DWs can be 
linearly stable, provided that the inter-site coupling in the lattice is sufficiently weak. 

From our results, it is clear that multi-component lattice models have a rich phenomenology, which is a natural 
addition to that of single-component ones. It would be interesting to observe the predicted features in experimental 
settings, including weakly coupled BECs and photonic-crystal nonlinear media. 
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